{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reactive Power Optimization with PandaModels.jl\n",
    "\n",
    "This tutorial describes how to solve some optimization problems regarding reactive power redispatch using **PandaModels**. The reactive power provision of generators, e.g., net.sgen, can be optimized for  \n",
    "\n",
    "1) maintaining voltage setpoints  \n",
    "2) maintaining reactive power setpoints  \n",
    "3) (active) power loss reduction \\\n",
    "4) branch loading reduction\n",
    "\n",
    "\n",
    "For more details on the corresponding models, see:\n",
    "* [maintaining voltage setpoints](https://github.com/e2nIEE/PandaModels.jl/blob/develop/src/models/vstab.jl)\n",
    "* [maintaining reactive power setpoints](https://github.com/e2nIEE/PandaModels.jl/blob/develop/src/models/qflex.jl)\n",
    "* [power loss reduction](https://github.com/e2nIEE/PandaModels.jl/blob/develop/src/models/ploss.jl)\n",
    "* [branch loading reduction](https://github.com/e2nIEE/PandaModels.jl/blob/develop/src/models/loading.jl)\n",
    "\n",
    "**Note**:\n",
    "\n",
    "- **PandaModels** is still in the development stage, and the developed models for reactive power optimizations are tested with limited grids. It is possible that errors may arise when facing different grids.\n",
    "- It is suggested to use **net.sn_mva=1.0** to start the optimization. Using other sn_mva values may introduce not converged OPF  \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1) Maintaining Voltage Setpoints\n",
    "\n",
    "Grid operators want to maintain a constant voltage profile across all or a part of the system buses. Reactive power optimization, in this case, helps them identify the preferable provision of reactive power from generators. The objective is to minimize the deviation of the voltages of certain buses (here, buses of DER) from a setpoint, e.g.,$V_{\\rm setpoint}\n",
    "= 1.03 p.u.$ As a result, the reactive power setpoint for each generator is calculated. The mathematical formulation of this objective function can have the following forms:\n",
    "\n",
    "\\begin{align}\n",
    "& \\underset{\\mathcal{X} = [q, ...]}{\\text{minimize}}\n",
    "& &  \\sum_{i\\in \\mathcal{BI}} (V_i-V_{\\rm setpoint})^2 \\\\\n",
    "& \\text{subject to}\n",
    "& & g(\\mathcal{X})=0 \\\\\n",
    "& & & h(\\mathcal{X}) \\leq 0\n",
    "\\end{align}\n",
    "\n",
    "where $V_{i}$ is the voltage variable of bus $i$ in $\\mathcal{BI}$ which denotes the set of buses of DER. The $g(\\mathcal{X})$ and $h(\\mathcal{X})$, denote equality and inequality constraints, respectively. The $\\mathcal{X}$ denotes the set of variables decisions, such as reactive power $q$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-1: Create test grid\n",
    "\n",
    "How does it work? Let's first create a cigre grid with DERs(distribution energy resource) from the pandapower's network database. After running the power flow calculation, we can see that there are nine DERs (net.sgen) in the grid, and they don't provide reactive power for normal operation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from copy import deepcopy\n",
    "from pandapower.networks.cigre_networks import create_cigre_network_mv\n",
    "from pandapower.run import runpp\n",
    "from pandapower.runpm import runpm_vstab, runpm_multi_vstab, runpm_qflex, runpm_ploss, runpm_loading"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "net = create_cigre_network_mv(with_der=\"pv_wind\")\n",
    "runpp(net) # run power flow calculation\n",
    "net_org = deepcopy(net) # copy of the grid for further comparison\n",
    "display(net.res_sgen) # show power flow results for DER"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(You can use the pandapower plotting module to visualize the grid model.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import pandapower.plotting as plot\n",
    "plot.simple_plot(net)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### 1-2: Parametrization for optimization\n",
    "Before starting the optimization, we need to configure some parameters. \n",
    "In PandaModels, the uncontrollable loads and DER (net.sgen) are converted to load, and the controllable power elements are converted to generators. Accordingly, we neet to set loads as uncontrollable and set DER (net.sgen) as controllable elements. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "net_opt = deepcopy(net) # copy of the grid for further comparison \n",
    "net_opt.load['controllable'] = False\n",
    "net_opt.sgen['controllable'] = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Now, let's set the constraints."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# lower and upper bounds for buses\n",
    "net_opt.bus[\"max_vm_pu\"] = 1.1\n",
    "net_opt.bus[\"min_vm_pu\"] = 0.9\n",
    "\n",
    "# lower and upper bounds for external grid\n",
    "net_opt.ext_grid[\"max_q_mvar\"] = 10000.0\n",
    "net_opt.ext_grid[\"min_q_mvar\"] = -10000.0\n",
    "net_opt.ext_grid[\"max_p_mw\"] = 10000.0\n",
    "net_opt.ext_grid[\"min_p_mw\"] = -10000.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "We keep the active power for DER constant to the calculated value from the power flow and define a constraint for the reactive power. (In this grid, there is no generator (net.gen))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# lower and upper bounds for DERs\n",
    "net_opt.sgen[\"max_p_mw\"] = net_opt.sgen.p_mw.values\n",
    "net_opt.sgen[\"min_p_mw\"] = net_opt.sgen.p_mw.values\n",
    "net_opt.sgen[\"max_q_mvar\"] = net_opt.sgen.p_mw.values * 0.328\n",
    "net_opt.sgen[\"min_q_mvar\"] = -net_opt.sgen.p_mw.values * 0.328\n",
    "display(net_opt.sgen)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Subsequently, let's set a high upper bound for lines and transformers to avoid congestion issues."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# lower and upper bounds for lines\n",
    "net_opt.trafo[\"max_loading_percent\"] = 100.0\n",
    "net_opt.line[\"max_loading_percent\"] = 100.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Finally, we need to let PandaModels know which buses are of interest and what is user-defined voltage setpoint. To do this, we can add extra column called **\"pm_param/setpoint_v\"** in net.bus and set values for buses connecting DER (net.sgen)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "net_opt.bus[\"pm_param/setpoint_v\"] = None # add extra column\n",
    "net_opt.bus.loc[net_opt.sgen.bus, \"pm_param/setpoint_v\"] = 0.96 # voltage setpoints for theses buses are 0.96"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-3: Start optimization\n",
    "Now let's call the function ***pandapower.runpm_vstab()*** to start the optimizaiton for maintaining voltage setpoints. \n",
    "(Note that the first time the function is called, Julia is started in the background, which may take some time.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "try:\n",
    "    runpm_vstab(net_opt)\n",
    "except Exception as err:\n",
    "    print(err)    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Also, the following parameters and options can be modified as inputs while calling the optimization from PandaModles. For this case, we can keep the default values.\n",
    "\n",
    "| parameter | description | type | default |\n",
    "| :--- | :--- | :---: | :--- |\n",
    "| correct_pm_network_data | checks if network data is correct. If not tries to correct it | bool | True |\n",
    "| silence | Suppresses information and warning messages output by PowerModels | bool | True |\n",
    "| pm_model | PowerModels.jl model to use | str | \"ACPPowerModel\" |\n",
    "| pm_solver | \"main\" solver| str | \"ipopt\" |\n",
    "| pm_mip_solver | mixed integer solver| str | \"cbc\" |\n",
    "| pm_nl_solver | nonlinear solver| str | \"ipopt\" |\n",
    "| pm_tol | default desired convergence tolerance for solver to use | float | 1e-8 |\n",
    "| pm_log_level | solver log level in power models | int | 0 |\n",
    "| delete_buffer_file | If True, the .json file used by PandaModels will be deleted after optimization. | bool | True |\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-4: Results comparison\n",
    "Let's check the results for voltage and reactive power at buses connecting DER before and after the optimization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# results from power flow\n",
    "display(net_org.res_bus.vm_pu[net_org.sgen.bus])\n",
    "v_org = net_org.res_bus.vm_pu[net_org.sgen.bus].values\n",
    "print(\"--- deviation from setpionts v=0.96 ---\")\n",
    "display(sum((v_org-0.96)**2)) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# results from optimization\n",
    "display(net_opt.res_bus.vm_pu[net.sgen.bus])\n",
    "v_opt = net_opt.res_bus.vm_pu[net_opt.sgen.bus].values\n",
    "print(\"--- deviation from setpionts v=0.96 ---\")\n",
    "display(sum((v_opt-0.96)**2)) # devition from setpoints v=0.96 p.u."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As exected, the bus voltages after optimization are close to the defined setpoints. It is caused by the provision of reactive power of DER."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# results from power flow\n",
    "display(net_org.res_sgen)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# results from optimization\n",
    "display(net_opt.res_sgen)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-5: Saving reactive power with multi-objectives\n",
    "In some situations, users want to improve the voltage profiles with limited reactive power provisions. For this, the objective function is modified by minimizing the reactive power provision. As the formulation below shows, the two sub-objectives for **maintaining voltage setpoints** and **saving reactive power provisions** are connected through coefficients $c_{\\rm v}$ and $c_{\\rm q}$. \n",
    "\n",
    "\\begin{align}\n",
    "& \\underset{\\mathcal{X} = [q, ...]}{\\text{minimize}}\n",
    "& &  \\sum_{i\\in \\mathcal{BI}} c_{\\rm v}*(V_i-V_{\\rm setpoint})^2 + c_{\\rm q}*(q_i-0)^2 \\\\\n",
    "\\end{align}\n",
    "\n",
    "We can call the objective function with different coefficient combinations, e.g., $c_{\\rm v}=0.99$ and $c_{\\rm q}=0.01$.\n",
    "(In default, $c_{\\rm v}=1.0$ and $c_{\\rm q}=0.0$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "net_optlim = deepcopy(net_opt)\n",
    "net_optlim[\"obj_factors\"] = [0.99, 0.01] # set coefficients for the multi-objective function\n",
    "try:\n",
    "    runpm_vstab(net_optlim)\n",
    "except:\n",
    "    print(\"[WinError 3] The system cannot find the path specified to the python39.dll\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As the results below show, with $c_{\\rm v}=0.99$ and $c_{\\rm q}=0.01$, the provided reactive powers are reduced, and the corresponding voltages are increased."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "print (\"--- optimized reactive power without limitation Cv = 1.0, Cq = 0. ----\")\n",
    "display(net_opt.res_sgen) # show optimized reactive power without limitation\n",
    "print (\"--- optimized reactive power with limitation Cv = 0.99, Cq = 0.01 ----\")\n",
    "display(net_optlim.res_sgen) # show optimized reactive power with limitation\n",
    "\n",
    "print (\"--- deviation from setpoints without limitation Cv = 1.0, Cq = 0. ----\")\n",
    "display(sum((v_opt - 0.96) ** 2)) # show optimized reactive power without limitation\n",
    "print (\"--- deviation from setpoints with limitation Cv = 0.99, Cq = 0.01 ----\")\n",
    "v_optlim = net_optlim.res_bus.vm_pu[net_optlim.sgen.bus].values\n",
    "display(sum((v_optlim - 0.96) ** 2)) # show optimized bus voltages"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The coefficient combination is grid-dependent. For example, the picture below shows the result for one of our case studies with a realistic grid. The user can find the most suitable combination setting through several attempts.\n",
    "<img src=\"pics/coef_combination.png\" alt=\"ALT\" width=\"600\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-6: Time series optimizations\n",
    "There are two ways to perform time series optimizations. \n",
    "The first one is optimization in a loop, i.e., updating the grid condition for each time step and repeatedly calling optimization. The second one uses the \"multi-timestep\" optimization. PowerModels can solve the multi-grid optimization problem. This feature is used in PandaModels to deal with time series optimization for on specific grid by replacing the \"multi-grid\" with \"multi-timestep\" of the grid. \n",
    "\n",
    "Next, you will see an example."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's load the time series data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "# load time series data for 96 time steps\n",
    "time_series = pd.read_json(\"cigre_timeseries_15min.json\")\n",
    "time_series.sort_index(inplace=True)\n",
    "\n",
    "# show time series\n",
    "display(time_series) \n",
    "\n",
    "# set load type and sgen type in the cigre grid\n",
    "net.load.loc[:, \"type\"] = \"residential\"\n",
    "net.sgen.loc[:, \"type\"] = \"pv\"\n",
    "net.sgen.loc[8, \"type\"] = \"wind\"\n",
    "\n",
    "# get rated power\n",
    "load_p = net.load.loc[:, \"p_mw\"].values\n",
    "sgen_p = net.sgen.loc[:7, \"p_mw\"].values\n",
    "wind_p = net.sgen.loc[8, \"p_mw\"]\n",
    "\n",
    "# create dataframe for loads and sgens\n",
    "load_ts = pd.DataFrame(index=time_series.index.tolist(), columns=net.load.index.tolist())\n",
    "sgen_ts = pd.DataFrame(index=time_series.index.tolist(), columns=net.sgen.index.tolist())\n",
    "for t in range(96):\n",
    "    load_ts.loc[t] = load_p * time_series.at[t, \"residential\"]\n",
    "    sgen_ts.iloc[t, 0:8] = sgen_p * time_series.at[t, \"pv\"]\n",
    "    sgen_ts.loc[t, 8] = wind_p * time_series.at[t, \"wind\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we perform the time series optimizations in a loop."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "########################################################################\n",
    "####################### Optimization in a loop #########################\n",
    "########################################################################\n",
    "\n",
    "y_loop = []\n",
    "for t in range(96):\n",
    "    \n",
    "    # Update operating points for each time step\n",
    "    net_opt.load.p_mw = load_ts.loc[t]\n",
    "    net_opt.sgen.p_mw = sgen_ts.loc[t]\n",
    "    \n",
    "    # If required, we also need to update the constraints, e.g., Qmax- and Qmin-value, for each time step. \n",
    "    net_opt.sgen.max_p_mw = sgen_ts.loc[t]   # for load: max_p_mw == min_p_mw, because they are not controllable.\n",
    "    net_opt.sgen.min_p_mw = sgen_ts.loc[t]\n",
    "    \n",
    "    # run optimization for time step t\n",
    "    try:\n",
    "        runpm_vstab(net_opt)\n",
    "    except:\n",
    "        print(\"[WinError 3] The system cannot find the path specified to the python39.dll\")\n",
    "    \n",
    "    # save the mean voltage value for each time step\n",
    "    y_loop.append(net_opt.res_bus.vm_pu[net_opt.sgen.bus].values.mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, we can do it using multi-timestep optimization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "########################################################################\n",
    "####################### multi-timestep optimization ####################\n",
    "########################################################################\n",
    "\n",
    "from pandapower.control import ConstControl\n",
    "from pandapower.timeseries import DFData\n",
    "\n",
    "# For this method, we need to add the time series data into the pandapower.controller.\n",
    "net_opt.controller.drop(net.controller.index, inplace=True) # remove the existing controllers\n",
    "\n",
    "# Add time series for load and sgen.\n",
    "ConstControl(net_opt, element=\"load\", variable=\"p_mw\",\n",
    "             element_index=net.load.index.tolist(), profile_name=net.load.index.tolist(),\n",
    "             data_source=DFData(load_ts))\n",
    "ConstControl(net_opt, element=\"sgen\", variable=\"p_mw\",\n",
    "             element_index=net.sgen.index.tolist(), profile_name=net.sgen.index.tolist(),\n",
    "             data_source=DFData(sgen_ts))\n",
    "\n",
    "# If required, we also need to creat a time series for the constraints, e.g., Qmax- and Qmin-value, for each time step. \n",
    "df_qmax, df_qmin = sgen_ts.copy(), sgen_ts.copy()\n",
    "df_qmax[df_qmax.columns] = net_opt.sgen.max_q_mvar\n",
    "df_qmin[df_qmin.columns] = net_opt.sgen.min_q_mvar\n",
    "\n",
    "# Then create controllers for them.\n",
    "ConstControl(net_opt, element=\"sgen\", variable=\"max_p_mw\",\n",
    "              element_index=net.sgen.index.tolist(), profile_name=net.sgen.index.tolist(),\n",
    "              data_source=DFData(sgen_ts))\n",
    "ConstControl(net_opt, element=\"sgen\", variable=\"min_p_mw\",\n",
    "              element_index=net.sgen.index.tolist(), profile_name=net.sgen.index.tolist(),\n",
    "              data_source=DFData(sgen_ts))\n",
    "ConstControl(net_opt, element=\"sgen\", variable=\"max_q_mvar\",\n",
    "              element_index=net.sgen.index.tolist(), profile_name=net.sgen.index.tolist(),\n",
    "              data_source=DFData(df_qmax))\n",
    "ConstControl(net_opt, element=\"sgen\", variable=\"min_q_mvar\",\n",
    "              element_index=net.sgen.index.tolist(), profile_name=net.sgen.index.tolist(),\n",
    "              data_source=DFData(df_qmin))\n",
    "\n",
    "# For the multi-timestep optimization, we need to call another PandaModels-function.\n",
    "try:\n",
    "    runpm_multi_vstab(net_opt, from_time_step=0, to_time_step=96)\n",
    "except:\n",
    "    print(\"[WinError 3] The system cannot find the path specified to the python39.dll\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note that**, currently, the function ***runpm_multi_vstab()*** doesn't support using of the multi-objective function for saving reactive power.\n",
    "\n",
    "The results for each time step are saved in **net.res_ts_opt**. Let's get them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "y_multi = []\n",
    "for t in range(96):\n",
    "    y_multi.append(net_opt.res_ts_opt[str(t)].res_bus.vm_pu[net_opt.sgen.bus].values.mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And then plot them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2, sharex=True, sharey=True)\n",
    "ax1.plot(range(96), y_loop)\n",
    "ax1.set_title(\"results for \\n optimization in a loop\")\n",
    "ax2.plot(range(96), y_multi)\n",
    "ax2.set_title(\"results for \\n multi-timestep optimization\")\n",
    "\n",
    "for ax in ax1, ax2:\n",
    "    ax.grid(True)\n",
    "    ax.label_outer()\n",
    "    ax.set_xlabel(\"time step [15min]\")\n",
    "ax1.set_ylabel(\"mean voltage value for DER-buses [p.u.]\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As the plot above shows, both methods have the same results.\n",
    "\n",
    "**For a long-term time series optimization (>500 timesteps), we suggest using the first method, optimization in a loop.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2) Maintaining Reactive Power Setpoints\n",
    "\n",
    "The motivation for this objective is to minimize the deviations of reactive power injection at the TSO-DSO interface from a target value, e.g., $Q_{\\rm setpoint}= 2 \\rm MVar$. It is formulated by the following equation:\n",
    "\n",
    "\\begin{align}\n",
    "& \\underset{\\mathcal{X} = [q, ...]}{\\text{minimize}}\n",
    "& &  \\sum_{i\\in \\mathcal{interfaces}} (Q_i-Q_{\\rm setpoint})^2 \\\\\n",
    "\\end{align}\n",
    "\n",
    "We use the same grid model to introduce the example. As the picture low shows, the marked bus (or interface) is observed.\n",
    "\n",
    "\n",
    "<img src=\"pics/qflex.png\" alt=\"ALT\" width=\"300\">\n",
    "\n",
    "First, let's create the grid and check the reactive power injection at this interface."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "net = create_cigre_network_mv(with_der=\"pv_wind\")\n",
    "runpp(net)\n",
    "net_org = deepcopy(net)\n",
    "\n",
    "display(net.res_trafo.q_lv_mvar[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Without optimization, the reactive power injection at this interface is -5.942 MVAr.\n",
    "\n",
    "Let's configure optimization parameters and define a target value for this interface, $Q_{\\rm setpoint}= -5 \\rm MVAr$, i.e., we optimize the reactive power provision of DER to achieve this target."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "net_opt = deepcopy(net) # copy of the grid for further comparison \n",
    "net_opt.load['controllable'] = False\n",
    "net_opt.sgen['controllable'] = True\n",
    "\n",
    "# lower and upper bounds for buses\n",
    "net_opt.bus[\"max_vm_pu\"] = 1.1\n",
    "net_opt.bus[\"min_vm_pu\"] = 0.9\n",
    "\n",
    "# lower and upper bounds for external grid\n",
    "net_opt.ext_grid[\"max_q_mvar\"] = 10000.0\n",
    "net_opt.ext_grid[\"min_q_mvar\"] = -10000.0\n",
    "net_opt.ext_grid[\"max_p_mw\"] = 10000.0\n",
    "net_opt.ext_grid[\"min_p_mw\"] = -10000.0\n",
    "\n",
    "# lower and upper bounds for DERs\n",
    "net_opt.sgen[\"max_p_mw\"] = net_opt.sgen.p_mw.values\n",
    "net_opt.sgen[\"min_p_mw\"] = net_opt.sgen.p_mw.values\n",
    "net_opt.sgen[\"max_q_mvar\"] = net_opt.sgen.p_mw.values * 0.328\n",
    "net_opt.sgen[\"min_q_mvar\"] = -net_opt.sgen.p_mw.values * 0.328\n",
    "display(net_opt.sgen)\n",
    "\n",
    "# lower and upper bounds for lines\n",
    "net_opt.trafo[\"max_loading_percent\"] = 100.0\n",
    "net_opt.line[\"max_loading_percent\"] = 100.0\n",
    "\n",
    "######################################################\n",
    "####################### Attention ####################\n",
    "######################################################\n",
    "\n",
    "# We need to let PandaModels know which interface is of interest and what is user-defined setpoint. \n",
    "# So, we need to add extra column called **\"pm_param/setpoint_q\"** in net.trafo and set target value. \n",
    "\n",
    "net_opt.trafo[\"pm_param/setpoint_q\"] = None # add extra column\n",
    "net_opt.trafo[\"pm_param/setpoint_q\"][0] = -5 #  reactive power setpoint for trafo 0 is -5 MVAR\n",
    "\n",
    "# Also, we need to tell PandaModels which side is observed.\n",
    "net_opt.trafo[\"pm_param/side\"] = None\n",
    "net_opt.trafo[\"pm_param/side\"][0] = \"lv\" # the observed side is the low-voltage side\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After parametrization, we can start the optimization using ***runpm_qflex()***."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "try:\n",
    "    runpm_qflex(net_opt)\n",
    "except:\n",
    "    print(\"[WinError 3] The system cannot find the path specified to the python39.dll\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's check the optimization results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "display(net_opt.res_trafo.q_lv_mvar[0])  # show Q-injection at the interface after optimization\n",
    "display(net_opt.res_sgen) # show optimized reactiv power provision of DER"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the results above, we can see that the Q-injection at the interface is closer to the setpoint after the optimization but not achieved. The reason is that the Q-maximum is limited. We can modify the Q-maximum and try again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "net_opt.sgen[\"max_q_mvar\"] = net_opt.sgen.p_mw.values * 0.628\n",
    "net_opt.sgen[\"min_q_mvar\"] = -net_opt.sgen.p_mw.values * 0.628\n",
    "try:\n",
    "    runpm_qflex(net_opt)\n",
    "except:\n",
    "    print(\"[WinError 3] The system cannot find the path specified to the python39.dll\")\n",
    "display(net_opt.res_trafo.q_lv_mvar[0])  # show Q-injection at the interface after optimization\n",
    "display(net_opt.res_sgen)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, the Q-injection at the interface is optimized to -5 MVAr. \n",
    "\n",
    "Similar to the **maintaining voltage setpoints**, if you are interested in time series optimizations, you can do it through a loop or use ***runpm_multi_qflex()***."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3) Power Loss Reduction and Branch Loading Reduction\n",
    "\n",
    "The corresponding objective functions are formulated as follows: \n",
    "\n",
    "- loss reduction\n",
    "\\begin{align}\n",
    "& \\underset{\\mathcal{X} = [q, ...]}{\\text{minimize}}\n",
    "& &  \\sum_{(i, j)\\in \\mathcal{N}} (S_{ij}+S_{ji})^2 \\\\\n",
    "\\end{align}\n",
    "\n",
    "\n",
    "- loading reduction\n",
    "\\begin{align}\n",
    "& \\underset{\\mathcal{X} = [q, ...]}{\\text{minimize}}\n",
    "& &  \\sum_{(i, j)\\in \\mathcal{N}} (i_{i->j, target}-0)^2 \\\\\n",
    "\\end{align}\n",
    "\n",
    "We still use the cigre medium voltage grid to introduce the case studies. First, let's create the grid, add constraints, and define the observed branches."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "net = create_cigre_network_mv(with_der=\"pv_wind\")\n",
    "net.load['controllable'] = False\n",
    "net.sgen['controllable'] = True\n",
    "net.sgen[\"max_p_mw\"] = net.sgen.p_mw.values\n",
    "net.sgen[\"min_p_mw\"] = net.sgen.p_mw.values\n",
    "net.sgen[\"max_q_mvar\"] = net.sgen.p_mw.values * 0.328\n",
    "net.sgen[\"min_q_mvar\"] = -net.sgen.p_mw.values * 0.328\n",
    "net.bus[\"max_vm_pu\"] = 1.1\n",
    "net.bus[\"min_vm_pu\"] = 0.9\n",
    "net.ext_grid[\"max_q_mvar\"] = 10000.0\n",
    "net.ext_grid[\"min_q_mvar\"] = -10000.0\n",
    "net.ext_grid[\"max_p_mw\"] = 10000.0\n",
    "net.ext_grid[\"min_p_mw\"] = -10000.0\n",
    "net.trafo[\"max_loading_percent\"] = 100.0\n",
    "net.line[\"max_loading_percent\"] = 100.0\n",
    "\n",
    "# add an additional column named \"pm_param/target_branch\" to let PandaModels know the observed branches.\n",
    "net.line[\"pm_param/target_branch\"] = True  \n",
    "net.switch.loc[:, \"closed\"] = True\n",
    "\n",
    "# run power flow calculation\n",
    "runpp(net)\n",
    "net_loss = deepcopy(net)\n",
    "net_loading = deepcopy(net)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you see, the power flow results without optimization are stored in the grid ***net***, and the grids ***net_loss*** and ***net_loading*** are going to used for optimization. Now we call the functions ***runpm_ploss*** and ***runpm_loading*** to reduce the power loss and line loading percent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "runpm_ploss(net_loss)\n",
    "runpm_loading(net_loading)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compare the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "print(\"loss without optimization:\", net.res_line.pl_mw.values.sum())\n",
    "print(\"loss with optimization:\", net_loss.res_line.pl_mw.values.sum())\n",
    "print(\"++++++++++++++++++++++++++++++++++++++++++++++\")\n",
    "print(\"loading without optimization:\", net.res_line.loading_percent.values.sum())\n",
    "print(\"loading with optimization:\", net_loading.res_line.loading_percent.values.sum())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also, we can reduce loss or loading percent of transformers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "net_loss.line.drop(columns=[\"pm_param/target_branch\"], inplace=True)\n",
    "net_loss.trafo[\"pm_param/target_branch\"] = True\n",
    "runpm_ploss(net_loss)\n",
    "print(\"loss without optimization:\", net.res_trafo.pl_mw.values.sum())\n",
    "print(\"loss with optimization:\", net_loss.res_trafo.pl_mw.values.sum())"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "execution": {
   "allow_errors": true,
   "timeout": 300
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
